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Abstract 

Three-dimensional  direct  numerical  simulation  (DNS)  of  exhaust  gas  recirculation  (EGR)-type  turbu¬ 
lent  combustion  operated  in  moderate  and  intense  low-oxygen  dilution  (MILD)  condition  has  been  carried 
out  to  study  the  flame  structure  and  flame  interaction.  In  order  to  achieve  adequate  EGR-type  initial/inlet 
mixture  fields,  partially  premixed  mixture  fields  which  are  correlated  with  the  turbulence  are  carefully  pre- 
processed.  The  chemical  kinetics  is  modelled  using  a  skeletal  mechanism  for  methane-air  combustion.  The 
results  suggest  that  the  flame  fronts  have  thin  flame  structure  and  the  direct  link  between  the  mean  reaction 
rate  and  scalar  dissipation  rate  remains  valid  in  the  EGR-type  combustion  with  MILD  condition.  How¬ 
ever,  the  commonly  used  canonical  flamelet  is  not  fully  representative  for  MILD  combustion.  During 
the  flame-flame  interactions,  the  heat  release  rate  increases  higher  than  the  maximum  laminar  flame  value, 
while  the  gradient  of  progress  variable  becomes  smaller  than  laminar  value.  It  is  also  proposed  that  the 
reaction  rate  and  the  scalar  gradient  can  be  used  as  a  marker  for  the  flame  interaction. 
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1.  Introduction 

Practical  combustion  systems  are  continuously 
required  to  be  more  efficient  and  more  environ¬ 
mentally  friendly.  In  conventional  combustion 
techniques,  preheating  the  unburnt  mixture  using 
exhaust  gases  is  one  way  to  improve  thermal  effi¬ 
ciency  of  the  system.  However,  preheating  also 
increases  the  flame  temperature,  resulting  in  an 
increase  of  thermal  NO  formation.  Therefore  it 
is  not  possible  to  achieve  both  objectives  using 
conventional  combustion  techniques. 

The  MILD  combustion  or  “flameless”  com¬ 
bustion  is  characterized  by  highly  preheated  mix¬ 
tures  and  low  temperature  rise  due  to 
combustion  [1-4],  in  which  (1)  reactants  are 
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diluted  with  large  amount  of  burnt  gases  so  that 
the  maximum  temperature  rise,  AT,  is  low 
compared  to  the  autoignition  temperature,  Tign, 
for  the  given  fuel  and  (2)  reactants  are  signifi¬ 
cantly  preheated,  higher  than  Tig„.  A  rigorous 
definition  of  MILD  regime  is  Tr  >  Tign  and 
AT  =  Tp  —  Tr  <  Tign  [4],  where  Tr  is  the  reactant 
temperature  and  Tp  is  the  product  temperature, 
and  is  shown  in  Fig.  1 .  Here,  the  autoignition  tem¬ 
perature  is  about  1100K  for  CH4-air  mixture 
(equivalence  ratio,  <j>  =  0.8),  which  is  calculated 
using  a  well  stirred  reactor  (WSR)  configuration 
with  a  residence  time  of  1  s.  This  definition  of 
MILD  combustion  is  based  on  WSR  theory  and 
is  therefore  easily  applied  to  premixed  combus¬ 
tion,  whereas  this  is  not  straightforward  for  non- 
premixed  systems. 

In  MILD  conditions,  the  flame  temperature  is 
very  low  compared  to  conventional  flames  due 
to  the  intense  dilution  even  with  the  highly 
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Fig.  1.  Combustion-type  diagram  [4]  showing  the 
condition  of  DNS  flame. 


preheated  unburnt  mixture.  For  instance.  Flame 
IV  shown  in  Fig.  1  has  a  preheat  temperature  of 
1 500  K,  while  its  flame  temperature  is  about 
1860  K.  Earlier  study  [1]  notes  that  a  sufficient  res¬ 
idence  time  for  thermal  NO  formation  is  some  sec¬ 
onds  at  around  1900  K  and  a  few  milliseconds  at 
around  2300  K,  which  clearly  suggests  that  the 
MILD  combustion  has  possibilities  for  significant 
NO*  reduction.  Therefore,  by  using  MILD  tech¬ 
nique,  it  is  possible  to  avoid  adverse  effect  of  pre¬ 
heating  on  flame  temperature  and  reduction  of 
NO*  emissions  is  no  longer  limited  by  the 
enhancement  of  thermal  efficiency. 

The  uniformity  of  temperature  due  to  the  high 
rate  of  recirculation  also  helps  to  reduce  combus¬ 
tion  instabilities  [1,2].  A  schematic  illustration  of 
combustion  with  EGR  is  shown  in  Fig.  2.  In  a 
combustor  employing  the  MILD  and  EGR  tech¬ 
niques,  a  portion  of  exhaust  gas  is  recirculated 
into  the  mixing  chamber  (denoted  as  ‘a’  in 
Fig.  2,  which  will  be  explained  in  Section  2.2  later) 
in  order  to  mix  with  fresh  reactants  and  exchange 
the  exhaust  heat  [5],  The  exhaust  and  fresh  gases 
may  be  partially  premixed  before  combustion 
occurs  in  the  combustion  chamber  (denoted  as 
‘b’  in  Fig.  2),  since  the  mixing  time  is  not  long 
enough  to  achieve  perfect  mixing  in  practice. 


Exhaust  gas 


Therefore,  combustion  that  takes  place  in  such 
systems  is  considered  to  be  different  from  the  tra¬ 
ditional  turbulent  flames:  exhaust  gas  pockets  of 
products  and  radicals,  which  are  not  well  mixed 
before  combustion,  can  lead  to  additional  com¬ 
plexities  such  as  flame-flame  interactions. 

Several  studies  have  been  carried  out  to  further 
our  understanding  of  MILD  combustion  [5-12], 
From  an  experimental  point  of  view,  no  flame 
front  is  visible  in  direct  photographs  of  MILD 
combustion  [6,9,12],  Plessing  et  al.  [5]  uses  laser 
diagnostics  to  compare  the  instantaneous  OH- 
PLIF  and  Rayleigh  thermometry  results  between 
conventional  and  MILD  combustion.  In  the 
MILD  case,  the  flame  front  seems  distributed 
from  the  temperature  field  images,  while  the  vari¬ 
ation  of  OH-PLIF  suggests  presence  of  flame 
fronts,  although  the  intensity  of  OH-PLIF  is  low 
compared  to  the  conventional  flame.  By  compari¬ 
son,  OH-PLIF  and  temperature  show  similar 
variation  in  conventional  combustion  conditions. 
Ozdemir  and  Peters  [9]  also  show  consistent 
results  from  comparison  between  temperature 
and  OH-PLIF  fields.  Conventional  numerical 
approaches  have  also  been  applied  to  MILD 
flames  [7,8,10,11]  using  either  flamelet  or  eddy- 
dissipation  based  methods  in  a  Reynolds  averaged 
Navier-Stokes  (RANS)  context.  In  most  RANS 
studies,  the  mean  velocity  and  temperature  fields 
show  a  consistent  trend  with  the  experimental 
data  obtained  by  laser  measurements.  However, 
the  predicted  peak  temperature  tends  to  differ 
from  the  experimental  results  in  the  range  of 
about  100-600  K,  especially  when  the  oxygen  con¬ 
centration  is  low  [11],  It  is  also  shown  by  Aminian 
et  al.  [7]  that  the  prediction  of  minor  species  such 
as  OH  and  CO  is  sensitive  to  the  temperature  fluc¬ 
tuation  in  RANS. 

Given  the  environmentally-friendly  nature  of 
this  combustion  mode,  it  is  very  useful  to  pose 
the  question:  what  is  the  flame  front  structure  in 
MILD  combustion?  We  believe  that  finding  an 
answer  to  this  question  would  help  to  construct 
a  modelling  framework  for  MILD  combustion. 
Specifically,  we  like  to  ask  whether  flamelet 
assumptions  are  valid  or  whether  the  flame  front 
is  still  flamelet  like  in  MILD  combustion  condi¬ 
tions:  from  a  direct  photograph  and  temperature 
measurements  shown  in  [5,6,9],  the  flamelet 
assumption  does  not  seem  to  be  valid;  whereas 
several  numerical  studies  [8,10]  using  flamelet 
models  show  comparable  results  with  experiments 
across  a  range  of  conditions. 

In  this  study,  a  three-dimensional  DNS  of 
CH4-air  partially  premixed  flame  under  MILD 
condition  is  carried  out  using  a  skeletal  mecha¬ 
nism  [13],  The  flame  is  analysed  from  several  view 
points  to  answer  the  questions  posed  above.  The 
DNS  detail  is  discussed  in  Section  2,  the  results 
are  presented  in  Section  3  and  the  conclusions 
are  drawn  in  Section  4. 
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2.  DNS  of  EGR-type  combustion 

2.1.  Configuration  and  numerical  implementation 

The  numerical  code  used  in  this  study  is 
SENGA2,  an  updated  version  of  SENGA  [14]. 
Compressible  transport  equations  are  solved  on 
a  uniform  grid  for  mass,  momentum,  total  inter¬ 
nal  energy,  and  the  mass  fraction  of  N  —  1  chem¬ 
ical  species  using  temperature  dependent 
transport  properties.  Elere,  a  skeletal  mechanism 
[13]  comprising  A=16  species  is  used.  Spatial 
derivatives  are  obtained  using  a  tenth  order  cen¬ 
tral  difference  scheme  which  gradually  reduces  to 
a  fourth  order  central  difference  scheme  near 
boundaries  and  fourth  order  one-sided  differenc¬ 
ing  on  boundaries.  Time  integration  is  achieved 
using  a  third  order  Runge-Kutta  scheme. 

Figure  3(a)  shows  the  numerical  configuration 
and  coordinates.  The  domain  is  cubic  with  a  non¬ 
reflecting  outflow  employing  Navier-Stokes  char¬ 
acteristic  boundary  conditions  (NSCBC)  [15]  on 
the  down  stream  (x-direction)  face,  and  periodic 
conditions  in  the  y  and  r-directions.  A  mixture 


Progress  variable 


Fig.  3.  (a)  Iso-surfaces  of  the  reaction  rate  and  x-y,  y-z 
and  z-x  planes  of  temperature  field.  The  grey,  transpar¬ 
ent  iso-surfaces  correspond  to  co+  =  0.7  and  coloured 
iso-surfaces  are  =  1.05,  in  which  interacting  and 
non-interacting  local  flame  fronts  are  respectively 
denoted  by  purple  and  light  blue  colours,  (b)  PDF  of 
progress  variable  at  x\,  x3,  x6,  x9  locations.  Blue  line: 
p(cY)  (initial  and  inlet  fields),  black  lines:  p(cT)  and  red 
lines:  p{cY).  (For  interpretation  of  the  references  to  color 
in  this  figure  legend,  the  reader  is  referred  to  the  web 
version  of  this  article.) 


of  exhaust  gas  and  fresh  premixed  gas  is  fed  from 
the  left  boundary  in  the  x-direction  at  an  average 
velocity  of  The  mass  fraction  of  partially  pre¬ 
mixed  mixture,  7„  and  turbulent  velocity,  h,  at  the 
inlet  ( x  =  0)  are  specified  as: 

u(x  =  0,y,z,t)  =  u{x{t),y,z),  (1) 

Yfx  =  0  ,y,z,t)  =  Yi(x{t),y,z),  (2) 

where  u  and  7,  correspond  to  the  velocity  and 
partially  premixed  scalar  fields  obtained  during 
the  preprocessing  stage  explained  in  Section  2.2. 
The  x-location  of  the  scanning  plane  at  time  t,  de¬ 
noted  by  x(r),  moves  with  velocity  through 
pre-computed  homogeneous  isotropic  turbulence 
and  scalar  fields,  and  it  and  7,  are  obtained  by 
interpolating  values  from  the  pre-computed  do¬ 
main  onto  the  scanning  plane. 

2.2.  Preprocessing  of  scalar  and  velocity  fields  for 
initial  and  inflow  conditions 

DNS  of  a  complete  EGR  combustion  system  is 
not  yet  feasible  because  of  heavy  computational 
cost.  Therefore,  the  combustion  phase  (b  in 
Fig.  2)  is  simulated  in  the  present  study,  and  the 
mixing  phase  (a  in  Fig.  2)  is  taken  into  account 
while  generating  the  initial  and  inflow  fields.  The 
mixture  used  in  this  DNS  is  partially  premixed 
between  fresh  reactants  and  the  exhaust  gases, 
which  is  assumed  as  the  inlet  mixture  for  the 
EGR-type  combustion.  The  steps  described  below 
are  followed  to  achieve  the  desired  fields  of  it  and 
7,  : 

(i)  The  turbulence  field  is  generated  in  a  preli¬ 
minary  DNS  of  freely  decaying,  homoge¬ 
neous  isotropic  turbulence  in  a  periodic 
domain  using  the  initial  turbulence  obtained 
as  in  [16],  and  the  simulation  is  continued 
until  the  turbulence  is  fully  developed. 

(ii)  A  one-dimensional  laminar  flame  is  simu¬ 
lated  in  a  desired  MILD  condition.  In  the 
ID  simulation,  the  unburnt  mixture  is 
diluted  with  H20  and  C02  in  such  way  that 
the  molar  proportion  of  H20  and  C02 
becomes  2:1  (which  is  the  case  for  complete 
combustion),  and  the  molar  concentration 
of  02  matches  the  desired  dilution  level. 

(iii)  A  turbulent  scalar  field  is  obtained  using 
scalar-energy  spectrum  function  as  in  [17]. 
This  field  is  taken  as  an  initial  progress  var¬ 
iable,  cr  field  with  values  between  0  and  1. 
Here,  cY=  1  —  YfiYfr,  where  Yf  denotes 
the  fuel  mass  fraction  and  the  subscript,  r, 
denotes  values  of  reactant  mixture.  In  order 
to  construct  a  mass  fraction  field  from  this 
cY  field,  the  ID  flame  result  obtained  in 
Step  (ii)  is  used.  The  mass  fraction  of  each 
species  is  tabulated  relative  to  the  progress 
variable  using  the  ID  flame  result,  and  the 
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temperature  is  set  to  a  constant  value  T'h,  to 
be  specified  later.  This  scalar  field  is  par¬ 
tially  premixed  between  reactants  and  prod¬ 
ucts,  and  the  scalar  fluctuations  do  not  yet 
have  any  correlation  with  the  turbulence 
obtained  in  Step  (i). 

(iv)  These  scalar  and  turbulence  fields  are  then 
allowed  to  evolve  in  a  periodic  domain  to 
mimic  the  EGR-mixing  without  reaction. 
The  duration  of  this  mixing  DNS  is  about 
one  large  eddy  timescale,  lju',  which  is 
much  shorter  than  the  autoignition  delay 
time,  where  /0  and  u'  are  respectively  the 
integral  length  scale  and  the  root- 
mean-square  (rms)  of  fluctuation  of  the  tur¬ 
bulence  obtained  in  Step  (i).  During  this 
procedure,  the  scalar  fields  develop  correla¬ 
tions  with  the  turbulence  field.  Temperature 
is  also  allowed  to  evolve  from  a  uniform 
value,  T'h  in  this  step  resulting  in  the  maxi¬ 
mum  temperature  fluctuation  of  about 
2.2%  of  the  mean  value.  The  mixture 
fraction  also  has  variations  within  a  range 
of  about  1.9%  of  the  mean  value,  which  is 
approximately  0.052  giving  tj)  =  0.94.  Here, 
the  mixture  fraction  is  calculated  using  the 
definition  proposed  in  [18], 

The  velocity  and  scalar  fields  obtained  in  the 
above  steps  are  used  as  the  initial  and  inlet  condi¬ 
tions  (a  in  Fig.  2)  in  the  combustion  DNS.  The 
probability  density  function  (PDF)  of  cY  in  the 
initial  field  (after  the  above  preprocessing)  is 
shown  in  Fig.  3(b)  (blue  line).  Although  following 
step  (iii),  the  progress  variable  field  has  a  typical 
bimodal  distribution,  with  two  sharp  peaks  at 
cY=  0  and  1,  during  step  (iv)  turbulent  mixing 
and  molecular  diffusion  produces  samples  with 
intermediate  values,  0  1,  with  significant 

PDFs  as  shown  in  Fig.  3(b).  The  mean  and  vari¬ 
ance  of  this  cr  field  are  respectively  (cY)  =  0.51 
and  (cy)  =  0.10. 

2.3.  Computational  parameters  and  conditions 

The  regime  of  combustion  in  the  present  DNS 
is  denoted  as  Flame  IV  in  Fig.  1.  The  equivalence 
ratio  of  the  injected  mixture  shown  in  Fig.  2, 
is  0.8  and  the  inlet  and  initial  mixture  tempera¬ 
tures  are  set  as  T'b  ~  1500  K.  This  inlet  tempera¬ 
ture  might  be  too  high  for  practical  systems. 
However  this  high  inlet  temperature  together  with 
the  dilution  level  used  in  this  study  shows  that  the 
flame  condition  is  strictly  in  the  MILD  regime  as 
one  can  see  in  Fig.  1 .  Moreover,  the  value  used  in 
this  study  is  comparable  to  that  used  by  Suzuka- 
wa  et  al.  [19],  The  maximum  molar  fraction  of 
oxygen  in  the  inlet  mixture,  X0lJ,  which  indicates 
the  dilution  level,  is  0.047,  and  the  averaged 
oxygen  molar  fraction,  (Xo2j),  is  0.035.  The 


maximum  flame  temperature  of  the  ID  laminar 
flame  used  in  step  (ii)  of  Section  2.2  is  about 
1860  K,  which  yields  AT  =  360  K.  The  unstrained 
laminar  flame  speed,  SL,  and  thermal  thickness, 
5th  =  (Tp  —  Tr)/\$T\max,  for  the  mixture  are 
respectively  3.20  m/s  and  0.691  mm.  The  Zeldo- 
vich  thickness,  SF,  which  is  defined  from  the  ratio 
between  the  thermal  diffusivity  and  laminar  flame 
speed  is  0.116  mm.  The  turbulence  field  obtained 
in  step  (iv)  has  u'/SL  =  5.13  and  l0j5F=  12.8, 
where  u'  is  the  turbulence  intensity  and  l0  is  the 
integral  length  scale.  The  turbulent  Reynolds 
numbers  based  on  l0  and  Taylor  length  scale,  2, 
are  Re/o  =  93.6  and  Re^  =  31.7  respectively.  The 
ratio  of  integral  length  scales  of  the  turbulence 
and  scalar  fields  is  about  1.58,  and  the  inlet  veloc¬ 
ity  Uin  =  7.825%  The  Karlovitz  number  is  esti¬ 
mated  using  Ka  (if j SL)3/2(l0/ ftf)112  and  the 
Damkohler  number  is  defined  as  Da  =  (10/SF)/ 
(if/SL).  In  the  present  DNS,  Ka  =  3.25  and 
Da  =  2.49,  which  is  classified  as  thin-reaction- 
zones  if  one  uses  the  classical  regime  diagram 
for  turbulent  premixed  flames  [20]. 

The  domain  has  dimensions  Lx  =  Lv  =  L-  = 
10.0  mm  (14.5<5(/l),  and  is  discretized  using 
512  x  512  x  512  uniform  grid  points,  ensuring  at 
least  20  grid  points  in  5th  (based  on  the  diagonal 
distance  between  grid  points). 

The  computational  domain  for  the  reacting 
flow  simulation  was  initialized  using  the  fields 
constructed  in  Section  2.2  and  the  simulation 
was  run  for  1.5  flow-through  times,  which  is  the 
mean  convection  time  from  the  inlet  to  outflow 
boundaries,  to  ensure  that  the  initial  transients 
had  left.  The  simulation  was  then  continued  for 
one  additional  flow-through  time  and  80  data  sets 
were  collected. 

DNS  has  been  conducted  on  Cray  XE6  system 
and  the  simulation  was  distributed  over  4096 
cores  and  256  nodes  requiring  a  wall  clock  time 
of  about  120  hours. 


3.  Results  and  discussion 

3.1.  DNS  results  and  flame  interactions 

Figure  3(a)  shows  iso-surfaces  of  the  reaction 
rate  and  x-y,  y-z  and  z-x  planes  of  temperature 
field.  The  grey,  transparent  iso-surfaces  of  reac¬ 
tion  rate  correspond  to  0%  =  0.7  and  coloured 
iso-surfaces  are  cnf  =  1.05,  in  which  interacting 
and  non-interacting  flame  elements  are  respec¬ 
tively  shown  as  purple  and  light  green  colours. 
The  criteria  to  distinguish  interacting  and 
non-interacting  flames  are  explained  later  in  this 
section.  The  superscript  +  denotes  appropriate 
normalization  using  pr,  SL  and  Sth,  where  pr  is 
the  reactant  mixture  density.  The  reaction  rate 
in  Fig.  3(a)  is  the  reaction  rate  of  a  temperature- 
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based  progress  variable,  ct=(T—  Tr)/{Tp—  Tr) 
and  so  is  equivalent  to  the  heat  release  rate.  One 
can  also  define  the  progress  variable  based  on 
c  Y-  Since  the  simulation  involves  non-unity  Lewis 
numbers,  these  two  progress  variables  will  be  dif¬ 
ferent  and  so  both  will  used  in  this  study. 

The  flame  shown  in  Fig.  3(a)  has  a  unique 
shape  compared  to  the  conventional  turbulent 
planar  flames  [21,22]  and  the  presence  of  flame 
interactions  is  obvious.  Due  to  the  preheat  tem¬ 
perature,  regions  of  intense  reaction  rate  are 
located  not  only  in  the  middle  of  computational 
domain  but  in  the  upstream  and  downstream 
regions  as  well.  The  maximum  temperature  in  this 
instantaneous  snap-shot  is  about  1710  K,  while  it 
is  about  1860  K  for  the  ID  laminar  flame  used  in 
the  step  (iv).  This  is  because  in  the  present  DNS, 
the  inlet  mixture  consists  of  both  the  diluted  pre¬ 
mixed  reactants  with  X0li  =  0.047  and  pockets 
of  their  products.  Thus  the  mean  oxygen  level 
becomes  lower  in  the  DNS  than  the  one-dimen¬ 
sional  laminar  calculation,  leading  to  lower  tem¬ 
perature  and  the  desirable  MILD  condition. 
Also,  in  conventional  turbulent  flames,  the  peak 
value  of  heat  release  rate  is  more  or  less  equal  to 
the  value  in  an  unstrained  laminar  flame 
(<u+  ~  1)  [21,22].  On  the  other  hand,  the  spatial 
extent  of  local  flames  having  co+  ~  1  in  the  present 
MILD  combustion  is  sparse  compared  to  conven¬ 
tional  turbulent  flames  [22],  as  shown  in  Fig.  3(a). 
This  is  caused  by  the  EGR  configuration,  and  is 
due  to  the  existence  of  exhaust  gas  pockets  and 
the  partially  premixed  state  of  inlet  mixture. 

Probability  density  function  of  the  progress 
variables,  cT  (black)  and  cY  (red),  is  shown  in 
Fig.  3(b)  for  four  x-locations:  x~[  =  1.25  (cT  = 
0.078,  cY  —  0.653);  =  3.88  (cj  =  0.188,  cY  = 
0.800);  x+  =  7.78  (cT  =  0.352,  cY  =  0.963); 
xj  =  11.7  ( cT  =  0.447,  cY  =  0.998).  Although  one 
can  distinguish  flame  fronts  in  terms  of  the  reac¬ 
tion  rate,  as  shown  in  Figs.  3(a)  and  Fig.  4,  the 
temperature  based  progress  variable  field  does 
not  have  a  clear  bimodal  distribution.  Especially, 


for  X3  location  the  PDF  of  cT  indicates  that  the 
probability  of  finding  reacting  gas  is  relatively 
high  (about  3-4).  This  is  because  the  flame  thick¬ 
ness  based  on  cT  in  the  DNS,  is  sometimes  much 
thicker  than  Sth  in  one-dimensional  DNS  due  to 
the  facts  that  (1)  heat  conduction  towards  the  ex¬ 
haust  gas  pockets  from  the  inlet,  which  do  not 
generate  much  heat  release  and  (2)  the  heat  release 
is  distributed  among  many  species,  and  so  is  more 
spread  out  in  space  because  of  differential  diffu¬ 
sion  compared  to  the  single  species  FCh4-  On  the 
other  hand,  the  PDFs  of  cjffred  lines)  show  differ¬ 
ent  distributions  compared  to  those  of  cT.  In  par¬ 
ticular,  the  probability  of  finding  fresh  reactants  is 
small  for  all  the  locations  shown  because  of  the 
presence  of  exhaust  gas  pockets  in  the  inlet  mix¬ 
ture,  and  the  mixing  and  fuel  consumption  in 
the  upstream  region.  However,  since  the  probabil¬ 
ity  of  finding  reacting  gas  is  also  small  compared 
to  the  PDF  of  cT,  it  is  observed  that  cY  field  has 
larger  gradient  and  so  thinner  fronts.  The  compar¬ 
ison  of  PDF  of  c T  and  cY  is  also  consistent  with 
the  difference  of  temperature  and  OH-PLIF  fields 
measured  in  past  experimental  studies  [5,9]  using 
EGR-type  combustor  in  which  the  OH-PLIF 
shows  clear  flame  fronts  while  the  flame  fronts 
in  the  temperature  field  seems  distributed.  If  one 
uses  unity  Lewis  number  in  DNS,  these  difference 
would  not  be  observed. 

Successive  snapshots  of  the  reaction  rate  are 
shown  in  Fig.  4(a)-(c)  in  a  x-z  plane  to  study 
the  flame  front  structure  and  flame  interactions. 
The  y-location  for  the  x-z  plane  changes  with  time 
depending  on  the  fluid  velocity  in  the  y-direction 
at  n+  =  0.  The  inset  figures  show  variations  of 
cT,  its  gradient  and  reaction  rate  as  a  function  of 
the  flame  normal  distance,  n+,  during  a  flame 
interaction  process  at  an  arbitrarily  chosen  loca¬ 
tion.  The  normal  distance  of  n+  =  0  at  each  time 
corresponds  to  the  circle  shown  on  the  reaction 
rate  contour  lines.  This  circle  moves  in  Lagrang- 
ian  sense  with  time.  It  should  be  noted  that  the 
normal  distance  is  calculated  using  the  local  \cT. 
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Fig.  4.  Two-dimensional  slice  of  reaction  rate  field.  Plots  of  an  arbitrarily  chosen  flame  interaction  process  as  a  function 
of  flame  normal  distance  is  shown  in  the  insets. 
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Therefore,  the  normal  direction  is  not  on  the  x-z 
plane  in  general.  The  duration  of  time  between 
Fig.  4(a)  and  4(c)  is  about  Q32if,  where  zF  = 
SJSL. 

It  is  generally  suggested  that  there  is  no  flame 
fronts  in  MILD  conbustion,  which  is  otherwise 
known  as  “flameless  combustion”  [1,4].  In  con¬ 
trast,  the  present  DNS  indicates  that  reactions 
occur  in  thin  flamelets  as  shown  in  Figs.  4(a)- 
(c).  However,  as  shown  in  Fig.  4(b),  top-left 
marked  with  a  dotted  square,  the  presence  of 
some  interacting  flames  gives  non-flamelet  like 
structure. 

The  inset  in  Fig.  4(a)  shows  that  two  local 
flames  have  a  thickness  of  about  5t/l  in  terms  of 
the  gradient  of  progress  variable.  The  flame  on 
the  left  has  a  flamelet  like  structure,  since  the  reac¬ 
tion  rate  is  sharply  peaked,  while  the  flame  on  the 
right  shows  a  relatively  large  thickness  and  smal¬ 
ler  gradients  of  cT.  The  inset  in  Fig.  4(b)  shows 
the  same  flames  as  in  the  inset  of  Fig.  4(a)  after 
0.16t^-  and  this  is  considered  as  the  beginning  of 
flame-flame  interaction.  At  this  stage  of  the  flame 
interaction,  the  reaction  rate  at  around  n+  =  0 
slightly  increases,  and  both  of  the  flames  have  a 
thin  flamelet  structure  with  peak  values  for  o>+ 
of  about  0.8  and  1.0.  Note  that  the  two  flame 
fronts  have  intense  heat  release  at  this  stage,  since 
cot  ~  1.  The  progress  variable  also  increases  at 
n+  =  0,  and  the  gradient  of  progress  variable 
shows  high  values  at  the  locations  of  intense  heat 
release.  This  is  typical  of  flamelet  type  combustion 
as  discussed  in  Section  3.2.  The  inset  in  Fig.  4(c) 
shows  that  two  flames  are  merged  in  terms  of 
the  reaction  rate.  The  peak  value  of  reaction  rate 
is  about  1.1,  which  is  larger  than  the  non-interact¬ 
ing  flames  shown  in  Fig.  4(b).  On  the  other  hand, 
the  normalized  gradient  of  progress  variable, 
|Vcr|+,  at  the  locations  of  intense  heat  release  rate 
(o>+  >  1  for  instance)  ranges  from  0.05  to  0.3, 
which  is  very  small  compared  to  the  non-interact¬ 
ing  flames  shown  in  the  inset  of  Fig.  4(b) 
(|Vct-|+  ~  0.5-0. 6  at  the  locations  of  intense  heat 
release  rate).  Thus  it  is  clear  that  conventional 
flamelet  type  approaches  might  not  be  able  to 
describe  the  flame-flame  interaction  sufficiently. 
Although,  the  maximum  value  of  the  normalized 
reaction  rate  during  the  flame  interaction  is  about 
1.1  in  Fig.  4,  values  up  to  about  1.6  (not  shown 
here)  were  found  in  the  data  in  interacting  flame 
regions.  Also,  the  locations  of  instantaneous  inter¬ 
acting  flames  are  shown  on  the  iso-surfaces  of 
reaction  rate  as  purple  colour  in  Fig.  3(a).  The 
criteria  for  flame  interaction  used  here  are 
(of  >  1.0  and  |Vc|+  <  0.3.  Since  the  reaction  rate 
of  the  interacting  flame  shows  non-interacting 
flame  like  structure,  it  is  not  possible  to  distin¬ 
guish  interacting  flames  using  only  the  reaction 
rate.  Although  these  thresholds  are  tentatively 
defined  based  on  Fig.  4(  a)— (c),  by  using  these  cri¬ 
teria  it  is  possible  to  distinguish  interacting  flames 


from  non-interacting  flames  in  order  to  assess  the 
effect  of  flame  interactions  on  models  in  future 
work. 

3.2.  Assessment  for  flamelet  approach 

It  is  also  our  interest  to  see  the  relation 
between  the  mean  reaction  rate,  a>c,  and  the  mean 
scalar  dissipation  rate,  ~tc,  especially  to  find  if 
flame  interactions  have  a  significant  effect  on  this 
relation  in  the  RANS  context.  Here  the  tilde 
denotes  Favre  (density-weighted)  time  average. 
Under  the  conditions  of  high  turbulent  Reynolds 
number  and  Damkohler  number  together  with 
unity  Lewis  number,  the  mean  reaction  rate  is 
related  to  the  mean  scalar  dissipation  rate  as  [23]: 

Wc  =  2C„,-1  ptc ’  (3) 

where  C„,  is  the  model  parameter  defined  as: 

Cm  —  [  ca>cP(c)dc/  f  cocP(c)dc,  (4) 

Jo  Jo 

and  the  typical  value  of  C,„  is  0.7-0. 8  for  lean  tur¬ 
bulent  premixed  hydrocarbon  flames,  where  P(c) 
is  the  marginal  PDF  of  the  progress  variable. 
Eq.  (3)  is  strictly  valid  when  the  flame  front  is  thin 
compared  to  the  Kolmogorov  length  scale  of  the 
turbulence.  Under  such  conditions,  the  progress 
variable  field  has  a  bimodal  distribution.  This  sit¬ 
uation  is  typically  known  as  flamelet  combustion 
in  general.  However,  this  expression  for  the  mean 
reaction  rate  is  sufficiently  accurate  even  for  thin 
reaction  zones  regime  combustion  [24],  In  order 
to  study  the  direct  relation  between  the  mean  reac¬ 
tion  rate  and  mean  scalar  dissipation  rate,  the 
model  parameter,  Cm,  is  calculated  from  Eq.  (3) 
for  both  cT  and  cY  and  is  shown  in  Fig.  5  as  a 
function  of  Favre-averaged  progress  variable.  As 
one  observes,  Cm  has  a  large  variation  for  cY, 
while  it  is  nearly  constant  for  cT ■  This  difference 
is  because  of  the  non-unity  Lewis  number  effects. 
Otherwise  these  two  values  for  Cm  would  be  close. 


Fig.  5.  Model  parameter  C,„  as  a  function  of  c. 
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However,  these  variations  of  C,„  supports  the  use 
of  Eq.  (3)  as  a  possible  model  for  the  mean  reac¬ 
tion  rate  in  MILD  combustion  also. 

By  denoting  the  volume  of  interacting  flame 
fronts  as  Fint,  obtained  using  the  criterion  given 
in  Section  3.1,  the  ratio  Vint/Vtot  is  about  36% 
when  the  flame  front  volume,  Vtot,  is  obtained 
using  the  condition  >  1.0.  This  ratio  drops 
to  about  5%  when  off  >  0.7  is  used  to  obtain  Vtot. 
Despite  these  observations,  the  relation  between 
the  mean  reaction  rate  and  the  scalar  dissipation 
rate  does  not  show  the  effect  of  the  flame  interac¬ 
tions  in  a  time-averaged  sense.  Also,  as  explained 
in  Section  3.1,  the  PDF  of  cT  shows  that  cT  field 
does  not  seem  to  have  a  bimodal  distribution. 
However,  it  can  be  said  based  on  Fig.  5  that  Eq. 
(3)  remains  sufficiently  valid.  These  results  are 
consistent  with  the  results  of  laser  measurements 

[5.9]  and  conclusions  of  RANS  simulations 

[8.10] ,  in  which  the  predicted  fields  have  consistent 
trends  with  experimental  results  while  measured 
temperature  fields  show  distributed  flame  fronts 
(non-bimodal  distribution  of  cT ). 

If  one  uses  Eq.  (4)  to  obtain  C„,  then  its  value 
would  strongly  depend  on  the  PDF  for  the  stan¬ 
dard  structure  of  the  flamelet  function  coc(c)  [23]. 
The  C,„  calculated  using  the  PDFs  in  Fig.  3(b) 
and  ®c(c)  obtained  from  the  unstrained  planar 
flamelets  is  shown  in  Fig.  5,  which  consists  the 
result  in  [23],  The  difference  in  C,„  values  obtained 
using  Eqs.  (3)  and  (4)  is  because  the  canonical 
flamelet  does  not  include  the  effects  of  flame  inter¬ 
actions  and  the  presence  of  exhaust  gas  pockets  in 
the  reactants.  This  suggests  that  an  appropriate 
canonical  flamelet  needs  to  be  identified  although 
the  flamelet  type  modelling  is  a  good  approxima¬ 
tion  for  the  MILD  combustion  condition  dis¬ 
cussed  here. 


4.  Conclusions 

Three-dimensional  DNS  of  EGR-type  turbu¬ 
lent  flame  has  been  conducted  under  a  MILD  con¬ 
dition  using  a  skeletal  mechanism  with  non-unity 
Lewis  numbers  and  following  conclusions  are 
obtained. 

•  The  presence  of  exhaust  gas  pockets  yields  a 
non-bimodal  PDF  for  cT  and  it  does  not  influ¬ 
ence  the  bimodal  behaviour  of  cY- 

•  The  reaction  are  observed  to  occur  in  thin 
regions  although  occasional  flame  interaction 
gives  a  distributed  structure. 

•  During  a  flame-flame  interaction,  the  local 
reaction  rate  becomes  larger  than  the 
unstrained  laminar  peak  value  while  the  gradi¬ 
ent  of  progress  variable  at  the  location  of 
intense  reaction  rate  becomes  very  small.  It  is 
also  proposed  that  these  two  quantities  can 
be  used  as  a  marker  for  the  flame  interactions. 


•  The  direct  relation  between  the  mean  reaction 
rate  and  scalar  dissipation  rate  is  observed  to 
be  sufficiently  valid  even  in  EGR-type  combus¬ 
tion  under  MILD  condition,  although  the  clas¬ 
sical  one-dimensional  unstrained  flame  might 
not  be  a  fully  representative  canonical  flamelet. 
The  construction  of  a  suitable  flamelet  is  of 
future  interest. 
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